#include "NeutrinoTools.h"

TString getTitleFromVar(const TString var)
{
  TString tit(var);
  tit.ReplaceAll("fMomErr/fMuonRec->P()","#sigma_{p}/p rec");
  tit.ReplaceAll("fNeutrinoSim->E()","E_{#nu sim} (GeV)");
  tit.ReplaceAll("fMuonSim->P()","p_{#mu sim} (GeV/c)");
  tit.ReplaceAll("fMuonSim->Pt()","p_{t}^{#mu sim int} (GeV/c)");
  tit.ReplaceAll("fMuonRec->Pt()","p_{t}^{#mu rec int} (GeV/c)");
  //tit.ReplaceAll("fMuonPtGuess->Mag()","p_{t}^{#mu rec, #nu sim} (GeV/c)");
  tit.ReplaceAll("fMuonSimPt->Mag()","p_{t}^{#mu sim} (GeV/c)");
  tit.ReplaceAll("fMuonRec->P()","p_{#mu rec} (GeV/c)");
  tit.ReplaceAll("fMuonRecPt->Mag()","p_{t}^{#mu rec} (GeV/c)");
  tit.ReplaceAll("fQ2Sim","Q^{2}_{sim} (GeV^{2}/c^{4})");
  tit.ReplaceAll("fMultiplicity","N_{all rec}");
  tit.ReplaceAll("fMuonRecFlightPath->Pt()","lever arm (m)");
  tit.ReplaceAll("fMuonRecFlightPath->Mag()","straight path length (m)");
  tit.ReplaceAll("(fMuonCharge[1]==-1)*2+(fMuonCharge[0]==-1)","TPC charge (11,-11,1-1,-1-1)");
  tit.ReplaceAll("fMuonRecNhits","selmu_tpc_nhits");
  tit.ReplaceAll("fMuonTypeSim","1:e, 2:#mu, 3:K_{L}, 4:#pi^{#pm}, 5:K^{#pm}, 6:p");
  tit.ReplaceAll("fNeutrinoParentType","1:e, 2:#mu, 3:K_{L}, 4:#pi^{#pm}, 5:K^{#pm}, 6:p");
  tit.ReplaceAll("fMuonAlphaRec","#theta_{#mu rec} (deg)");
  tit.ReplaceAll("fMuonAlphaSim","#theta_{#mu sim} (deg)");
  tit.ReplaceAll("fNeutrinoRec->Theta()*TMath::RadToDeg()","#theta_{#nu rec} (deg)");
  tit.ReplaceAll("fNeutrinoSim->Theta()*TMath::RadToDeg()","#theta_{#nu sim} (deg)");
  tit.ReplaceAll("fMuonRecEndPos->X()","endposX (m)");
  tit.ReplaceAll("fMuonRecEndPos->Y()","endposY (m)");
  tit.ReplaceAll("fMuonRecEndPos->Z()","endposZ (m)");
  tit.ReplaceAll("fMuonRecVertex->X()","vertexX (m)");
  tit.ReplaceAll("fMuonRecVertex->Y()","vertexY (m)");
  tit.ReplaceAll("fMuonRecVertex->Z()","vertexZ (m)");
  tit.ReplaceAll("fMuonRec->Theta()*TMath::RadToDeg()","#theta_{#mu rec lab} (deg)");
  tit.ReplaceAll("(fNeutMode==1||fNeutMode==2)","CCQE");
  tit.ReplaceAll("fMuonRecChi2/fMuonRecNDOF","#chi^{2}_{track}/NDOF");
}

plotMom(){
  TString tag="mom";
  gStyle->SetOptStat(0);

  t0=new TH1D("t0","",100,0.1,20); NeutrinoTools::BinLog(t0->GetXaxis()); t0->SetTitle("red: no cut, blue: M>1, green: M>1 & mu_P > 1 GeV/c ;Q^{2}_{sim} (GeV^2);Counts"); t0->SetLineColor(kRed); t0->Sumw2(); t0->SetLineWidth(2);
  t1=new TH1D("t1","",100,0.1,20); NeutrinoTools::BinLog(t1->GetXaxis()); t1->SetLineColor(kBlue); t1->Sumw2();t1->SetLineWidth(2);
  t2=new TH1D("t2","",100,0.1,20); NeutrinoTools::BinLog(t2->GetXaxis()); t2->SetLineColor(kGreen); t2->Sumw2();t2->SetLineWidth(2);
  
  tree->Draw("fQ2Sim>>t0");
  tree->Draw("fQ2Sim>>t1","fMultiplicity>1");
  tree->Draw("fQ2Sim>>t2","fMultiplicity>1 && fMuonRec.P()>1");

  e10=(TH1D*)t1->Clone("e10"); e10->Divide(t1,t0,1,1,"B"); e10->SetTitle("blue: M>1, green: M>1 & mu_P > 1 GeV/c, black: green/blue; Q^{2}_{sim} (GeV^2); #epsilon"); 
  e20=(TH1D*)t2->Clone("e20"); e20->Divide(t2,t0,1,1,"B");
  e21=(TH1D*)t2->Clone("e21"); e21->Divide(t2,t1,1,1,"B"); e21->SetLineColor(kBlack);
  
  c1=new TCanvas;
  gPad->SetLogx();
  gPad->SetLogy();
  t0->Draw();
  t1->Draw("same");
  t2->Draw("same");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.png");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.eps");

  c2=new TCanvas;
  gPad->SetLogx();
  e10->Draw("e");
  e20->Draw("same e");
  e21->Draw("same e");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.png");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.eps");
}

plotPt(){
  TString tag="pt";
  gStyle->SetOptStat(0);

  t0=new TH1D("t0","",100,0.1,20); NeutrinoTools::BinLog(t0->GetXaxis()); t0->SetTitle("red: no cut, blue: M>1, green: M>1 & mu_pT > 1 GeV/c ;Q^{2}_{sim} (GeV^2);Counts"); t0->SetLineColor(kRed); t0->Sumw2();t0->SetLineWidth(2);
  t1=new TH1D("t1","",100,0.1,20); NeutrinoTools::BinLog(t1->GetXaxis()); t1->SetLineColor(kBlue); t1->Sumw2();t1->SetLineWidth(2);
  t2=new TH1D("t2","",100,0.1,20); NeutrinoTools::BinLog(t2->GetXaxis()); t2->SetLineColor(kGreen); t2->Sumw2();t2->SetLineWidth(2);

  tree->Draw("fQ2Sim>>t0");
  tree->Draw("fQ2Sim>>t1","fMultiplicity>1");
  tree->Draw("fQ2Sim>>t2","fMultiplicity>1 && fMuonRecPt.Mag()>1");

  e10=(TH1D*)t1->Clone("e10"); e10->Divide(t1,t0,1,1,"B"); e10->SetTitle("blue: M>1, green: M>1 & mu_pT > 1 GeV/c, black: green/blue; Q^{2}_{sim} (GeV^2); #epsilon");
  e20=(TH1D*)t2->Clone("e20"); e20->Divide(t2,t0,1,1,"B");
  e21=(TH1D*)t2->Clone("e21"); e21->Divide(t2,t1,1,1,"B"); e21->SetLineColor(kBlack);

  c1=new TCanvas;
  gPad->SetLogx();
  gPad->SetLogy();
  t0->Draw();
  t1->Draw("same");
  t2->Draw("same");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.png");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.eps");

  c2=new TCanvas;
  gPad->SetLogx();
  e10->Draw("e");
  e20->Draw("same e");
  e21->Draw("same e");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.png");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.eps");
}

plotPtGuess(){
  TString tag="PtGuess";
  gStyle->SetOptStat(0);

  t0=new TH1D("t0","",100,0.1,20); NeutrinoTools::BinLog(t0->GetXaxis()); t0->SetTitle("red: no cut, blue: M>1, green: M>1 & mu_PtGuess > 1 GeV/c ;Q^{2}_{sim} (GeV^2);Counts"); t0->SetLineColor(kRed); t0->Sumw2();t0->SetLineWidth(2);
  t1=new TH1D("t1","",100,0.1,20); NeutrinoTools::BinLog(t1->GetXaxis()); t1->SetLineColor(kBlue); t1->Sumw2();t1->SetLineWidth(2);
  t2=new TH1D("t2","",100,0.1,20); NeutrinoTools::BinLog(t2->GetXaxis()); t2->SetLineColor(kGreen); t2->Sumw2();t2->SetLineWidth(2);

  tree->Draw("fQ2Sim>>t0");
  tree->Draw("fQ2Sim>>t1","fMultiplicity>1");
  tree->Draw("fQ2Sim>>t2","fMultiplicity>1 && fMuonPtGuess>1");

  e10=(TH1D*)t1->Clone("e10"); e10->Divide(t1,t0,1,1,"B"); e10->SetTitle("blue: M>1, green: M>1 & mu_PtGuess > 1 GeV/c, black: green/blue; Q^{2}_{sim} (GeV^2); #epsilon");
  e20=(TH1D*)t2->Clone("e20"); e20->Divide(t2,t0,1,1,"B");
  e21=(TH1D*)t2->Clone("e21"); e21->Divide(t2,t1,1,1,"B"); e21->SetLineColor(kBlack);

  c1=new TCanvas;
  gPad->SetLogx();
  gPad->SetLogy();
  t0->Draw();
  t1->Draw("same");
  t2->Draw("same");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.png");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.eps");

  c2=new TCanvas;
  gPad->SetLogx();
  e10->Draw("e");
  e20->Draw("same e");
  e21->Draw("same e");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.png");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.eps");
}


plotPtInv(){
  TString tag="ptinv";

  gStyle->SetOptStat(0);

  t0=new TH1D("t0","",100,0.1,20); NeutrinoTools::BinLog(t0->GetXaxis()); t0->SetTitle("red: no cut, blue: mu_pT>1, green: M>1 & mu_pT > 1 GeV/c ;Q^{2}_{sim} (GeV^2);Counts"); t0->SetLineColor(kRed); t0->Sumw2();t0->SetLineWidth(2);
  t1=new TH1D("t1","",100,0.1,20); NeutrinoTools::BinLog(t1->GetXaxis()); t1->SetLineColor(kBlue); t1->Sumw2();t1->SetLineWidth(2);
  t2=new TH1D("t2","",100,0.1,20); NeutrinoTools::BinLog(t2->GetXaxis()); t2->SetLineColor(kGreen); t2->Sumw2();t2->SetLineWidth(2);

  tree->Draw("fQ2Sim>>t0");
  tree->Draw("fQ2Sim>>t1","fMuonRecPt.Mag()>1");
  tree->Draw("fQ2Sim>>t2","fMultiplicity>1 && fMuonRecPt.Mag()>1");

  e10=(TH1D*)t1->Clone("e10"); e10->Divide(t1,t0,1,1,"B"); e10->SetTitle("blue: mu_pT>1, green: M>1 & mu_pT > 1 GeV/c, black: green/blue; Q^{2}_{sim} (GeV^2); #epsilon");
  e20=(TH1D*)t2->Clone("e20"); e20->Divide(t2,t0,1,1,"B");
  e21=(TH1D*)t2->Clone("e21"); e21->Divide(t2,t1,1,1,"B"); e21->SetLineColor(kBlack);

  c1=new TCanvas;
  gPad->SetLogx();
  gPad->SetLogy();
  t0->Draw();
  t1->Draw("same");
  t2->Draw("same");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.png");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.eps");

  c2=new TCanvas;
  gPad->SetLogx();
  e10->Draw("e");
  e20->Draw("same e");
  e21->Draw("same e");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.png");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.eps");
}

plotPtGuessInv()
{
  TString tag="PtGuessinv";

  gStyle->SetOptStat(0);

  t0=new TH1D("t0","",100,0.1,20); NeutrinoTools::BinLog(t0->GetXaxis()); t0->SetTitle("red: no cut, blue: mu_PtGuess>1, green: M>1 & mu_PtGuess > 1 GeV/c ;Q^{2}_{sim} (GeV^2);Counts"); t0->SetLineColor(kRed); t0->Sumw2();t0->SetLineWidth(2);
  t1=new TH1D("t1","",100,0.1,20); NeutrinoTools::BinLog(t1->GetXaxis()); t1->SetLineColor(kBlue); t1->Sumw2();t1->SetLineWidth(2);
  t2=new TH1D("t2","",100,0.1,20); NeutrinoTools::BinLog(t2->GetXaxis()); t2->SetLineColor(kGreen); t2->Sumw2();t2->SetLineWidth(2);

  tree->Draw("fQ2Sim>>t0");
  tree->Draw("fQ2Sim>>t1","fMuonPtGuess>1");
  tree->Draw("fQ2Sim>>t2","fMultiplicity>1 && fMuonPtGuess>1");

  e10=(TH1D*)t1->Clone("e10"); e10->Divide(t1,t0,1,1,"B"); e10->SetTitle("blue: mu_PtGuess>1, green: M>1 & mu_PtGuess > 1 GeV/c, black: green/blue; Q^{2}_{sim} (GeV^2); #epsilon");
  e20=(TH1D*)t2->Clone("e20"); e20->Divide(t2,t0,1,1,"B");
  e21=(TH1D*)t2->Clone("e21"); e21->Divide(t2,t1,1,1,"B"); e21->SetLineColor(kBlack);

  c1=new TCanvas;
  gPad->SetLogx();
  gPad->SetLogy();
  t0->Draw();
  t1->Draw("same");
  t2->Draw("same");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.png");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.eps");

  c2=new TCanvas;
  gPad->SetLogx();
  e10->Draw("e");
  e20->Draw("same e");
  e21->Draw("same e");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.png");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.eps");
}

/*

.L plot.C
double tol=0.05; TString tag = "muonRecPtNeutrinoE"; Double_t x0 = 0.09, x1 = 11; Bool_t klogx=1; while(tol<=0.21){ plotFlux(tol, tag, x0, x1, klogx); tol+=0.05;}
double tol=0.05; TString tag = "AlphaRecNeutrinoE"; Double_t x0 = 0, x1 = 180; Bool_t klogx=0; while(tol<=0.21){ plotFlux(tol, tag, x0, x1, klogx); tol+=0.05;}
double tol=0.05; TString tag = "muonPRecNeutrinoE"; Double_t x0 = 0.09, x1 = 15; Bool_t klogx=1; while(tol<=0.21){ plotFlux(tol, tag, x0, x1, klogx); tol+=0.05;}
double tol=0.05; TString tag = "MultiplicityNeutrinoE"; Double_t x0 = 0, x1 = 10; Bool_t klogx=0; while(tol<=0.21){ plotFlux(tol, tag, x0, x1, klogx); tol+=0.05;}

*/

plotFlux(const Double_t tol = 0.05, const TString tag="muonRecPtNeutrinoE", const Double_t x0=0.09, const Double_t x1=11, const Bool_t klogx=1)
{
  gStyle->SetOptStat(0);

  hh=(TH2D*) gDirectory->Get(tag);
  if(!hh){
    printf("no hist %s\n", tag.Data()); return;
  }

  TH2D *hmap=0x0;
  TH2D *hdiff=0x0;

  TGraphAsymmErrors *gr = NeutrinoTools::GetFluxMap(tol, hh, hmap, hdiff);
  
  c1=new TCanvas; gPad->SetLogx(); gPad->SetLogy(klogx);  gPad->SetLogz(); hh->Draw("colz"); c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"original.png");
  c1=new TCanvas; gPad->SetLogx(); gPad->SetLogy(klogx);  hmap->SetMinimum(0); hmap->SetMaximum(2); hmap->Draw("colz"); c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"map1-2.png");
  c1=new TCanvas; gPad->SetLogx(); gPad->SetLogy(klogx);  hmap->SetMinimum(0.5); hmap->SetMaximum(1.5); hmap->Draw("colz"); c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"map0.5-1.5.png");
  
  c1=new TCanvas; gPad->SetLogx(); gPad->SetLogy(klogx);  hdiff->SetMinimum(0); hdiff->SetMaximum(1); hdiff->Draw("colz"); c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"diff0-1.png");
  c1=new TCanvas; gPad->SetLogx(); gPad->SetLogy(klogx);  hdiff->SetMinimum(0); hdiff->SetMaximum(0.1); hdiff->Draw("colz"); c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"diff0-0.1.png");
  
  c1=new TCanvas; gPad->SetGrid();  gPad->SetLogx(klogx); gPad->SetLogy(); gr->SetMarkerSize(1); gr->SetMarkerStyle(24); gr->GetXaxis()->SetLimits(x0, x1); gr->GetYaxis()->SetRangeUser(0.2, 21); gr->Draw("alp"); c1->Print(Form("outplot/%scurveTol%.2f.png", tag.Data(), tol));

}

plotEMode(const Bool_t khigh=kTRUE)
{
  TCanvas *c1=new TCanvas;

  gStyle->SetOptStat(0);

  const Int_t mm = 30;
  const TString cut0 = "fMuonNSim>0 &&";
  const TString ecut = khigh? "fNeutrinoSim->E()>5.318&&fNeutrinoSim->E()<11.774" : "fNeutrinoSim->E()>0.490&&fNeutrinoSim->E()<1.085";

  TH1D * hh1 = new TH1D("hh1","",mm,-0.5, mm-0.5); hh1->SetLineWidth(4); 
  TH1D * hh2 = new TH1D("hh2","",mm,-0.5, mm-0.5); hh2->SetLineWidth(4); 

  //mode<0 for anti-neutrino
  tree->Draw("abs(fNeutMode)>>hh1",cut0+ecut+" && fQ2Sim<0.9"); hh1->SetLineColor(kRed);  hh1->SetTitle(ecut+", red: Q^{2}<0.9 GeV^{2}, blue: Q^{2}>1 GeV^{2}; Neut Mode; Counts");  
  tree->Draw("abs(fNeutMode)>>hh2",cut0+ecut+" && fQ2Sim>1");   hh2->SetLineColor(kBlue); 

  printf("hh1 ou %.0f hh2 ou %.0f\n", hh1->GetBinContent(0)+hh1->GetBinContent(hh1->GetNbinsX()), hh2->GetBinContent(0)+hh2->GetBinContent(hh2->GetNbinsX()));

  const Double_t hm = khigh? hh2->GetMaximum() : hh1->GetMaximum();

  hh1->SetMaximum( hm*1.1); hh1->SetMinimum(0);hh1->Draw();hh2->Draw("same");
  c1->Print(Form("outplot/neutmodelogy%dkhigh%d.png", gPad->GetLogy(), khigh));
  c1->Print(Form("outplot/neutmodelogy%dkhigh%d.eps", gPad->GetLogy(), khigh));
  
  gPad->SetLogy(1);
  hh1->SetMaximum( hm*2); hh1->SetMinimum(0.5);hh1->Draw();hh2->Draw("same");
  c1->Print(Form("outplot/neutmodelogy%dkhigh%d.png", gPad->GetLogy(), khigh));
  c1->Print(Form("outplot/neutmodelogy%dkhigh%d.eps", gPad->GetLogy(), khigh));
}

plotptMode()
{
  TCanvas *c1=new TCanvas;

  gStyle->SetOptStat(0);

  const Int_t mm = 30;
  TH1D * hh1 = new TH1D("hh1","",mm,-0.5, mm-0.5); hh1->SetLineWidth(5);
  TH1D * hh2 = new TH1D("hh2","",mm,-0.5, mm-0.5); hh2->SetLineWidth(5);
  TH1D * hh3 = new TH1D("hh3","",mm,-0.5, mm-0.5); hh3->SetLineWidth(4);
  TH1D * hh4 = new TH1D("hh4","",mm,-0.5, mm-0.5); hh4->SetLineWidth(3);

  tree->Draw("abs(fNeutMode)>>hh1","fMuonRecPt->Mag()<1"); hh1->SetLineColor(kRed);  
  tree->Draw("abs(fNeutMode)>>hh2","fMuonRecPt->Mag()>1"); hh2->SetLineColor(kBlue); 
  tree->Draw("abs(fNeutMode)>>hh3","fMuonRecPt->Mag()>1.5"); hh3->SetLineColor(kGreen); 
  tree->Draw("abs(fNeutMode)>>hh4","fMuonRecPt->Mag()>2"); hh4->SetLineColor(kMagenta); 
  const Double_t sf3 = hh2->GetBinContent(27)/hh3->GetBinContent(27);
  hh3->Scale(sf3);

  const Double_t sf4 = hh2->GetBinContent(27)/hh4->GetBinContent(27);
  hh4->Scale(sf4);

  const TString tit=Form("red: p_{t}^{#mu, rec}<1 GeV/c, blue: p_{t}^{#mu, rec}>1 GeV/c, green: #times %.1f p_{t}^{#mu, rec}>1.5 GeV/c, magenta: #times %.1f p_{t}^{#mu, rec}>2 GeV/c; Neut Mode; Counts", sf3, sf4);  
  hh1->SetTitle(tit);hh2->SetTitle(tit);

  printf("OU: hh1 %.0f hh2 %.0f hh3 %.0f hh4 %0.f\n", hh1->GetBinContent(0)+hh1->GetBinContent(hh1->GetNbinsX()), hh2->GetBinContent(0)+hh2->GetBinContent(hh2->GetNbinsX()), hh3->GetBinContent(0)+hh3->GetBinContent(hh3->GetNbinsX()), hh4->GetBinContent(0)+hh4->GetBinContent(hh4->GetNbinsX()));

  const Double_t hm = hh1->GetMaximum();

  hh2->Draw(); hh3->Draw("same"); hh4->Draw("same");
  c1->Print(Form("outplot/PtModelogy%dk.png", gPad->GetLogy()));
  c1->Print(Form("outplot/PtModelogy%dk.eps", gPad->GetLogy()));

  gPad->SetLogy(1);
  hh1->SetMaximum( hm*2); hh1->SetMinimum(0.5);hh1->Draw();hh2->Draw("same"); hh3->Draw("same"); hh4->Draw("same");
  c1->Print(Form("outplot/PtModelogy%dk.png", gPad->GetLogy()));
  c1->Print(Form("outplot/PtModelogy%dk.eps", gPad->GetLogy()));


}


plotNuE(const TString tag)
{
  ll=lout;

  TDirectory *tmpdir = new TDirectory("tmpdir","tmpdir");
  for(Int_t ii=0; ii<ll->GetEntries(); ii++){
    htmp = ll->At(ii);
    //printf("loading %d: %s\n", ii, htmp->GetName());
    tmpdir->Add(htmp);
  }

  tmpdir->cd();

  c1=new TCanvas;
  gStyle->SetOptStat(0);
  NeutrinoE->SetTitle((tag.Contains("/")?tag:tag.Prepend("outplot/"))+";E_{#nu, sim}; Counts");
  gPad->SetLogx();
  NeutrinoE->Draw();

  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"nuE.png");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"nuE.eps");

  delete tmpdir;
}

plotResPtGuess()
{
  gROOT->LoadMacro("style.cxx+");

  style::SetGlobalStyle();

  style::fgkYTitleOffset = 1.4;
  style::fgkXTitleOffset = 1.4;

  TCanvas *c1= new TCanvas;
  style::PadSetup(c1);

  Double_t lm=0.15;
  Double_t rm=0.13;
  Double_t tm=0.03;
  Double_t bm=0.17;

  c1->SetLeftMargin(lm);
  c1->SetRightMargin(rm);
  c1->SetTopMargin(tm);
  c1->SetBottomMargin(bm);

  TH2D * hh = new TH2D ("hh","",50,0.1, 2.5,50,-1,1); NeutrinoTools::BinLog(hh->GetXaxis());
  hh->SetTitle(";p_{t}^{#mu sim} (GeV/c);p_{t}^{#mu rec}/p_{t}^{#mu sim}-1");
  style::ResetStyle(hh);
  tree->Draw("fMuonPtGuess->Mag()/fMuonSimPt->Mag()-1:fMuonSimPt->Mag()>>hh","fMuonNSim>0 && fMuonNRec>0 &&  fMultiplicity>=3");
  TH2D *hnor = NeutrinoTools::NormalHist(hh, 5, 1);
  style::ResetStyle(hnor);
  gPad->SetLogx();

  hh->Draw("colz");
  //c1->Print("outplot/PtGuessRes.eps");
  c1->Print("outplot/PtGuessRes.png");

  hnor->Draw("colz");
  //c1->Print("outplot/PtGuessNorRes.eps");
  c1->Print("outplot/PtGuessNorRes.png");
}

//====================================================================================
//====================================================================================


void ptSpectrum(const Bool_t kqacut, const Int_t imode=0, const Bool_t kproton = kFALSE)
{
  const TString pwd=gSystem->pwd();

  if(pwd.Contains("SK")){
    if(kqacut){
      return;
    }
  }

  //===============================================

  gStyle->SetOptStat(0);
  
  gROOT->LoadMacro("style.cxx+");

  style::SetGlobalStyle();

  style::fgkYTitleOffset = 1.4;
  style::fgkXTitleOffset = 1.4;

  TCanvas *c0= new TCanvas;
  style::PadSetup(c0);

  gPad->SetLogy();
  gPad->SetLogx();

  Double_t lm=0.15;
  Double_t rm=0.02;
  Double_t tm=0.04;
  Double_t bm=0.17;

  c0->SetLeftMargin(lm);
  c0->SetRightMargin(rm);
  c0->SetTopMargin(tm);
  c0->SetBottomMargin(bm);

  gStyle->SetOptTitle(1);

  const Int_t nxbins = 15;
  const Double_t xmin = 0.05;
  const Double_t xmax = 2;

  //===============================================

  //========>
  TString xRecVar = kqacut?"fMuonRecPt->Mag()":"fMuonRec->Pt()";

  TString modename;
  const TString modecut = GetModeCut(imode, modename);
  const TString finalcut = Form("(fMuonNRec > 0 %s && %s)", kqacut? Form(" && %s",(GetQACut()).Data()): "", modecut.Data());
  if(kproton){
    xRecVar.ReplaceAll("Muon","Proton");
  }
  printf("xRecVar %s || finalcut { %s } || kqacut %d \n",xRecVar.Data(), finalcut.Data(), kqacut);
  //<=========

  TH1D * hrec = new TH1D("hrec","",nxbins, xmin, xmax); style::ResetStyle(hrec);   hrec->SetTitle(Form("%s;p_{t}^{#mu} (GeV/c);Counts",""));   
  NeutrinoTools::BinLog(hrec->GetXaxis()); 
  tree->Draw(xRecVar+">>hrec",finalcut);

  TH1D * hsim = new TH1D("hsim","",nxbins, xmin, xmax); style::ResetStyle(hsim); 
  NeutrinoTools::BinLog(hsim->GetXaxis()); 
  tree->Draw("fMuonSimPt->Mag()>>hsim",finalcut);

  TLegend * lg = new TLegend(0.17,0.8,0.25,0.95);
  style::ResetStyle(lg,0.3);
  lg->SetTextAlign(12);
  //lg->SetBorderSize(1);

  lg->SetHeader(Form("%s %s", kqacut?"D-V w/ QA":"#vec{#nu}=(001) w/o QA", modename.Data()));
  hrec->SetMinimum(TMath::Max(1.0,hsim->GetMinimum()*0.5));
  hrec->SetMaximum(hrec->GetBinContent(hrec->GetMaximumBin())*2);
  hrec->SetLineWidth(2);
  hrec->SetLineColor(kBlack);
  hrec->Draw("hist e"); 
  lg->AddEntry(hrec,"rec","l");

  hsim->SetLineWidth(3);
  hsim->SetLineColor(kRed);
  hsim->SetLineStyle(kDashed);
  hsim->Draw("hist e  same");
  lg->AddEntry(hsim,"sim","l");

  lg->Draw();
  TString tag=Form("kqacut%d_mode%s",kqacut, modename.Data());
  tag.ReplaceAll("#","");
  tag.ReplaceAll("{","");
  tag.ReplaceAll("}","");
  tag.ReplaceAll(" or ","");

  c0->Print(Form("outplot/ptSpec_%s.png",tag.Data()));


  c1->SetTopMargin(tm);
  gPad->SetLogx();
  gPad->SetLogy(0);

  gPad->SetGrid();
 
  hrec->Sumw2();
  hsim->Sumw2();
  hrec->Divide(hsim);
  hrec->SetYTitle("spectrum ratio (rec/sim)");
  hrec->SetMarkerStyle(20);
  hrec->SetMarkerSize(1);
  hrec->SetMinimum(0.5);
  hrec->SetMaximum(1.5);
  hrec->Draw("e");

  TLatex *lt=new TLatex(0.55,0.9, lg->GetHeader());
  style::ResetStyle(lt);
  //lt->SetTextAlign(33);
  lt->Draw();

  TF1 * f1=new TF1("f1","1",0,10);
  style::ResetStyle(f1);
  f1->SetLineColor(kRed);
  f1->SetLineWidth(2);
  f1->Draw("same");
  c0->Print(Form("outplot/ptSpRatio_%s.png",tag.Data()));

}

TF1* tmpplotFuncblock(TH1D *hh, TF1 *ftmp, const Int_t kcol, const Int_t lsty,  const TString fitopt, TLegend *lg)
{
  printf("\nfitting %s %s...\n", ftmp->GetName(), ftmp->GetTitle());

  style::ResetStyle(ftmp); 
  ftmp->SetParameters(100,0,0.1); 
  ftmp->SetLineColor(kcol); 
  ftmp->SetLineWidth(3);
  ftmp->SetLineStyle(lsty);
  hh->Fit(ftmp,fitopt); 
  lg->AddEntry(ftmp,ftmp->GetName(),"l");

  return ftmp;
}

plotFunc(const TString fin)
{
  new TFile(fin);

  TH2D *hori = h2;

  const Int_t ibin = 13;
  TH1D *hh = h2->ProjectionY("proj",ibin, ibin);
  const TString header(Form("%s %.2f - %.2f", h2->GetXaxis()->GetTitle(), h2->GetXaxis()->GetBinLowEdge(ibin), h2->GetXaxis()->GetBinUpEdge(ibin)));
  hh->SetTitle(h2_RMSmpv->GetTitle());
  hh->SetYTitle("Counts");
  hh->SetXTitle(hori->GetYaxis()->GetTitle());

  //===

  gROOT->LoadMacro("style.cxx+");

  style::SetGlobalStyle();

  style::fgkYTitleOffset = 1.4;
  style::fgkXTitleOffset = 1.4;

  TCanvas *c0= new TCanvas;
  style::PadSetup(c0);

  Double_t lm=0.15;
  Double_t rm=0.02;
  Double_t tm=0.13;
  Double_t bm=0.17;

  c0->SetLeftMargin(lm);
  c0->SetRightMargin(rm);
  c0->SetTopMargin(tm);
  c0->SetBottomMargin(bm);

  //===

  gStyle->SetOptTitle(1);
  gStyle->SetTitleW(0.99);
  style::ResetStyle(hh);
  hh->SetLineColor(kBlack);
  hh->SetMarkerColor(kBlack);
  hh->SetMarkerStyle(20);
  hh->SetMarkerSize(0.7);
  //gPad->SetLogy();
  hh->SetMaximum(hh->GetBinContent(hh->GetMaximumBin())*(gPad->GetLogy()?3:1.3));

  TLegend * lg = new TLegend(0.17,0.44,0.85,0.83);
  style::ResetStyle(lg,0.07);
  lg->SetTextAlign(12);
  //lg->SetBorderSize(1);
  
  //fits are affected by SetRangeUser!!
  TF1 * fgausML =   tmpplotFuncblock(hh, new TF1("GausML","[0]*TMath::Gaus(x,[1],[2],1)",-2,2),       kRed,     kDashed,     "NL", lg);
  TF1 * fgausLS =   tmpplotFuncblock(hh, new TF1("GausLS","[0]*TMath::Gaus(x,[1],[2],1)",-2,2),       kBlue,    kDotted,     "N",  lg);
  TF1 * fcauchyML = tmpplotFuncblock(hh, new TF1("CauchyML","[0]*TMath::CauchyDist(x,[1],[2])",-2,2), kGreen,   kSolid,      "NL", lg);
  TF1 * fcauchyLS = tmpplotFuncblock(hh, new TF1("CauchyLS","[0]*TMath::CauchyDist(x,[1],[2])",-2,2), kMagenta, kDashDotted, "N",  lg);

  //GetRMS is affected by SetRangeUser
  lg->SetHeader(Form("#splitline{RMS=%.2f, #sigma_{G ML}=%.2f, #sigma_{G LS}=%.2f, #sigma_{C ML}=%.2f, #sigma_{C LS}=%.2f}{%s}", hh->GetRMS(),  fgausML->GetParameter(2), fgausLS->GetParameter(2), fcauchyML->GetParameter(2), fcauchyLS->GetParameter(2),header.Data()));

  //======================================================

  //SetRangeUser only after getting all numbers
  hh->GetXaxis()->SetRangeUser(-0.65,0.65);

  hh->Draw("e");

  TLine * ln = new TLine(0,0,0,hh->GetMaximum());
  ln->SetLineColor(kGray);
  ln->SetLineWidth(1);
  ln->Draw("same");

  fgausML->Draw("same");
  fgausLS->Draw("same");
  fcauchyML->Draw("same");
  fcauchyLS->Draw("same");
  hh->Draw("same");



  lg->Draw();

  TString outn(fin);
  outn=outn(outn.First("/")+1,outn.Length());
  outn=outn(0,outn.First("."));
  cout<<outn<<endl;
  c0->Print(Form("outplot/momDescription%s%s.png", outn.Data(), hh->GetName()));

  //======================================================

  gPad->SetLogy(0);
  gPad->SetLogx();

  TH1D *hchi = 0x0;
  const TString hns[]={"GausML", "CauchyML", "GausLS", "CauchyLS"};
  for(Int_t ii=0; ii<4; ii++){
    TH1D * htmp=(TH1D*)gDirectory->Get(Form("h2_%schi", hns[ii].Data()));
    if(htmp){
      if(hchi){
        printf("\nbad!! hchi already exist! %s\n", hchi->GetName()); exit(1);
      }
      else{
        hchi=htmp;
      }
    }
  }
  if(!hchi){
    printf("\nbad!! null hchi!!\n"); exit(1);
  }

  hchi->SetMaximum(7.5);
  hchi->SetMinimum(0);
  hchi->Draw();

  TF1 * f1=new TF1("f1","1",0,10);
  style::ResetStyle(f1);
  f1->SetLineColor(kRed);
  f1->SetLineWidth(2);
  f1->Draw("same");

  c0->Print(Form("outplot/momDescription%schi.png", outn.Data()));
}

void subFit(TTree * tt,  TH2D *hh, const TString yVar, const TString xVar, const TString cut, const Int_t varopt, const Double_t fitmin, const Double_t fitmax)
{
  const TString tag(hh->GetName());

  TCanvas *cc= new TCanvas;
  style::PadSetup(cc);

  const Double_t lm=0.15;
  Double_t rm=0.14;
  const Double_t tm=0.13;
  const Double_t bm=0.17;

  cc->SetLeftMargin(lm);
  cc->SetRightMargin(rm);
  cc->SetTopMargin(tm);
  cc->SetBottomMargin(bm);

  gStyle->SetOptTitle(1);
  gStyle->SetOptFit(1);
  gStyle->SetTitleW(0.95);
  gStyle->SetTitleX(0.3);
  gStyle->SetStatW(0.15);
  gStyle->SetStatX(1);

  //===========================

  TList * lout=new TList;
  lout->Add(hh);

  //======
  const TString plvar(Form("%s:%s>>%s",yVar.Data(), xVar.Data(), hh->GetName()));

  printf("\n\ncheck plvar %s || cut { %s }\n", plvar.Data(), cut.Data());
  tt->Draw(plvar,cut);

  hh->SetTitle(Form(";%s;%s",getTitleFromVar(xVar).Data(), getTitleFromVar(yVar).Data()));
  hh->Draw("colz");
  cc->Print(Form("outplot/%s2d.png", tag.Data()));

  //============================================================
  rm = 0.03;
  cc->SetRightMargin(rm);

  const Double_t thres = 20;
  TH1D *hrmsnor=0x0, *hrmsmpv=0x0, *hrmswid=0x0, *hrmsres=0x0, *hrmschi=0x0;
  NeutrinoTools::FitSlicesY(hh, hrmsnor, hrmsmpv, hrmswid, hrmsres, hrmschi, "RMS", thres, lout);

  TH1D *hnor=0x0, *hmpv=0x0, *hwid=0x0, *hres=0x0, *hchi=0x0;
  NeutrinoTools::FitSlicesY(hh, hnor, hmpv, hwid, hres, hchi, "CauchyML", thres, lout);

  hwid->SetMarkerStyle(20);
  hwid->SetMarkerSize(1);
  hwid->SetMarkerColor(kBlack);
  hwid->SetLineWidth(1);
  hwid->SetLineColor(kBlack);

  hrmswid->SetMarkerStyle(22);
  hrmswid->SetMarkerSize(1.5);
  hrmswid->SetMarkerColor(kRed);
  hrmswid->SetLineColor(kRed);
  hrmswid->SetLineWidth(1);

  TH1D * hfit = 0x0;
  TH1D * hrms = 0x0;
  if(varopt==0){
    hfit=hrmsmpv;
    hfit->SetMaximum(1.2);
  }
  else{
    hfit=hwid;
    hfit->SetMaximum(0.6);
    //hrms=hrmswid;
    //hrms->SetMaximum(0.6);
  }
  
  const Double_t hmax = hrms? hrms->GetBinContent(hrms->GetMaximumBin()) : hfit->GetBinContent(hfit->GetMaximumBin());
  //hfit->SetMaximum(hmax*1.3);
  hfit->SetMinimum(0);
  TF1 * ff=new TF1("ff","[0]/pow(x,[1])", fitmin, fitmax);
  ff->SetParameters(0.2,1);
  if(hrms){
    hrms->Fit(ff, "R");
  }
  else{
    hfit->Fit(ff, "R");
  }
  hfit->SetTitle(Form("%s  %s",hh->GetName(), ff->GetTitle()));
  lout->Add(ff);

  style::ResetStyle(hfit);
 

  if(hrms){
    style::ResetStyle(hrms);
    hrms->Draw();
  }

  hfit->Draw(hrms?"same":"");

  cc->Print(Form("outplot/%svaropt%d.png", tag.Data(), varopt));

  //========

  TFile * fout=new TFile(Form("outplot/fitslice%s.root", tag.Data()),"recreate");
  lout->Write();
  fout->Save();
  fout->Close();
  delete fout;
}

void fitLeverArm(const TString fin, const Int_t mode)
{
  TFile::Open(fin);
  TTree *tt = tree;

  gROOT->LoadMacro("style.cxx+");

  style::SetGlobalStyle();

  style::fgkYTitleOffset = 1.4;
  style::fgkXTitleOffset = 1.4;

  //===================================================
  const TString basecut = "fMuonNRec>0 && fMuonNSim>0 && fMuonRec->P()<3.5";

  if(mode < 10){
    const Double_t nhitmin=55+mode*5;
    const Double_t nhitmax = nhitmin+5;
    const TString tag = Form("nhit%.0f%.0f", nhitmin, nhitmax);
    const TString nhitcut = basecut+Form(" &&  fMuonRecNhits[0]> %.0f && fMuonRecNhits[0]<%.0f", nhitmin, nhitmax);
    //subFit(tt, new TH2D(tag+"dd","", 12,0, 4, 60, -2, 2),  "(1/fMuonRec->P()-1/fMuonSim->P())/(1/fMuonRec->Pt()-1/fMuonSim->Pt())",   "fMuonRecFlightPath->Pt()", nhitcut, 1, 0.7, 3.1);
    subFit(tt, new TH2D(tag+"PToverP","", 12,0, 4, 110, 0, 1.1), "fMuonSim->Pt()/fMuonSim->P()",      "fMuonRecFlightPath->Pt()", nhitcut, 0, 0.7, 3.1);
    subFit(tt, new TH2D(tag+"dpT","", 12,0, 4, 60, -1.2, 1.2),   "1/fMuonRec->Pt()-1/fMuonSim->Pt()", "fMuonRecFlightPath->Pt()", nhitcut, 1, 0.7, 3.1);
    subFit(tt, new TH2D(tag+"dMom","", 12,0, 4, 60, -1.2, 1.2),  "1/fMuonRec->P()-1/fMuonSim->P()",   "fMuonRecFlightPath->Pt()", nhitcut, 1, 0.7, 3.1);
  }
  else{
    mode -=10;
    const Double_t larmmin=0.5+mode*0.5;
    const Double_t larmmax = larmmin+0.5;
    const TString tag = Form("larm%03.0f%03.0f", larmmin*100, larmmax*100);
    const TString larmcut = basecut+Form(" &&  fMuonRecFlightPath->Pt()> %.2f && fMuonRecFlightPath->Pt()<%.2f", larmmin, larmmax);
    
    subFit(tt, new TH2D(tag+"PToverP","", 10, 0, 100, 110, 0, 1.1), "fMuonSim->Pt()/fMuonSim->P()",      "fMuonRecNhits[0]", larmcut, 0, 30, 80);
    subFit(tt, new TH2D(tag+"dpT","", 10, 0, 100, 60, -1.2, 1.2),   "1/fMuonRec->Pt()-1/fMuonSim->Pt()", "fMuonRecNhits[0]", larmcut, 1, 30, 80);
    subFit(tt, new TH2D(tag+"dMom","", 10, 0, 100, 60, -1.2, 1.2),  "1/fMuonRec->P()-1/fMuonSim->P()",   "fMuonRecNhits[0]", larmcut, 1, 30, 80);
  }
}

cmpMomMuon()
{
  gStyle->SetOptStat(0);

  c1=new TCanvas;

  gPad->SetLogx();
  gPad->SetLogz();

  hh=new TH2D("hh","",100,0.05,2,300,0,3);NeutrinoTools::BinLog(hh->GetXaxis());
  
  tree->Draw("fTmpMom/fMuonRec->P():fMuonRecPt->Mag()>>hh","fMuonNRec>0 && ( fMuonRec->P()<3.5 && fMuonRecNhits[0]>60 && fMuonRecFlightPath->Pt()>1.5  ) && ( fMuonRecNhits[0]<80 && fMuonRecFlightPath->Mag()>3.0 )","colz"); 
  h2=NeutrinoTools::NormalHist(hh,2,1); h2->Draw("colz");
  h2->SetTitle("mom'/mom:pT, mom'=reconstructed with muon hypo. after QA cuts");
  h2->GetYaxis()->SetRangeUser(0.5,1.5);
  h2->Draw("colz");
  c1->Print("outplot/momMuon1.png");

  c1=new TCanvas;
  gPad->SetLogx();
  gPad->SetLogz();

  tree->Draw("fTmpMom/fMuonRec->P():fMuonRecPt->Mag()>>hh","fMuonNRec>0","colz"); 
  h2=NeutrinoTools::NormalHist(hh,2,1); h2->Draw("colz");
  h2->SetTitle("mom'/mom:pT, mom'=reconstructed with muon hypo. no QA cuts");
  h2->GetYaxis()->SetRangeUser(0.5,1.5);
  h2->Draw("colz");
  c1->Print("outplot/momMuon2.png");
}
